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Two elastically coupled nanomechanical resonators driven independently near their resonance 
frequencies show intricate nonlinear dynamics. The dynamics provide a scheme for realizing a 
nanomechanical system with tunable frequency and nonlinear properties. For large vibration ampli- 
tudes the system develops spontaneous oscillations of amplitude modulation that also show period 
doubling transitions and chaos. The complex nonlinear dynamics are quantitatively predicted by a 
simple theoretical model. 

PACS numbers: 05.45.-a, 62.25.-g, 85.85.+j 



Resonant nanoelectromechanical systems (NEMS) [ij 
are attracting interest in a broad variety of research areas 
and for many possible applications due to their remark- 
able combination of properties: small mass, high operat- 
ing frequency, large quality factor, and easily accessible 
nonlinearity [3| ■ Further development of NEMS applica- 
tions such as superior mass force f3| and charge [H| 
sensors, or reaching the quantum limit of detection in me- 
chanical systems [6|, requires addressing several impor- 
tant challenges. For example, the nonlinearity of the de- 
vices has to be either minimized or utilized for improving 
performance 0|- In addition, the large-scale integration 
of nanodevices demands a detailed understanding of the 
behavior of coupled devices in NEMS arrays [8, 9, 10, 11]. 

In this paper we demonstrate complex nonlinear be- 
havior of a pair of coupled nanomechanical devices, and 
show that this can be quantitatively understood from 
the basic physics of the devices. We show that the linear 
and weakly nonlinear response of one oscillation can be 
modified by driving the second oscillation, and, for some 
ranges of parameters of the devices, that the linear re- 
sponse range of the first oscillation can be significantly 
extended. When both oscillations are driven into their 
strongly nonlinear range more complicated frequency- 
sweep response curves are found, corresponding to the 
well known bistability of driven anharmonic "Duffing" 
resonators, but now with switching between a variety 
of different stable states of the coupled pair. Sponta- 
neous amplitude modulation oscillations may develop, 
with frequencies characteristic of the dissipation rates 
rather than of the intrinsic frequencies or their sums and 
differences. These amplitude modulations show period 
doubling bifurcations and chaos. The complex dynam- 
ics are reproduced quantitatively by a simple theoretical 
model, giving us confidence that the nonlinear behav- 
ior of coupled nanomechanical devices can be understood 
and controlled. 

We study a system of two strongly coupled nonlin- 
ear nanoelectromechanical resonators using a structure 
of doubly-clamped beams with a shared mechanical ledge 
shown in Fig. la. The devices consist of a stack of three 



layers of gallium arsenide (GaAs): a lOOnm highly n- 
doped layer, a 50nm insulating layer, and another 50nm 
layer that is highly p-dopcd. The piezoelectric property 
of GaAs results in a highly efficient integrated actuation 
mechanism described in 12]. A preliminary 120nm deep 
etch step is done to isolate the actuation electrodes of the 
two beams, so that the two beams can be addressed sep- 
arately while retaining strong elastic coupling. Optical 
interferometry is used for the motion transduction [l3| . 
The laser beam is adjusted so that both beams are in the 
illuminated spot. 

We model the behavior of the two strongly interacting 
nonlinear resonators by a system of coupled equations of 
motion for the beam displacements Xi,X2 in their funda- 
mental modes 

xi -t- 7iii -I- ujIxi -t- aix\ + D{xi — X2) — .gm(i), (la) 
X2 + 12X2 + ^1x2 + a2x\ + D{x2 - xi) = gD2{t). (lb) 

As well as the usual terms describing the resonant fre- 
quencies, damping, and Duffing nonlinearity, we include 
a linear coupling term in the displacements of strength 
D. The terms on the right hand side are the external 
drives applied to the two beams, which are controlled in- 
dependently. The linear terms in the equations, ignoring 
for now the drive and dissipation, give two modes with 
frequencies and ujn and corresponding eigenvectors 
e 1,^1 1 [HI- The frequency difference of the modes re- 
sults from both the intrinsic frequency difference ui — L02 
and the coupling D. 

To investigate the behavior of the system, the two 
beams are connected to two different sources with in- 
dependent frequencies and amplitudes. The monitored 
output variables are the amplitudes and phases of a lin- 
ear combination of the mechanical displacements of the 
two beams at or near the two drive frequencies f2]|. By 
applying various combinations of small amplitude sig- 
nals to the two beams at frequencies near the mode res- 
onances, the linear coupling parameters can be deter- 
mined. For the particular device shown in Fig. la the 
mode frequencies are determined to be oji/in — 16.79 
MHz and ujj/2tt = 17.25 MHz and the eigenvectors 
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FIG. 1: (color online) (a) SEM image of the system (beam 
dimensions: 6/im x 500nm x 200nm). (b) Up (dark, blue) and 
down (light, yellow) frequency sweeps of the amplitude |A| of 
the response for a single drive of strength 4.3 times critical 
as a function of the frequency relative to the linear resonance 
frequency luq and scaled by the width 7. The amplitude is 
plotted in units of the maximum amplitude at the critical 
drive strength, (c) Frequency shift of the weakly driven first 
mode as a function of the displacement amplitude \Aii \ of the 
strongly driven second mode: points and solid line - experi- 
ment; dashed line - small amplitude theory. 



e/ = (0.854,0.521) and e// = (-0.521,0.854). The fre- 
quencies of the individual resonators determined from in- 
verting the mode equations are uji/2tt — 16.71 MHz and 
u!2/2tt — 16.92 MHz. The frequency separation of 200 
kHz is consistent with the fabrication tolerance. The 
coupling strength \Id jln — 2.63 MHz is consistent with 
the strength of the elastic coupling found by finite ele- 
ment simulations of similar devices flBj. Transforming 
the measured width of the modes back to the original 
equations determines the values of 71 , 72 . The nonlinear 
parameters a\^ai are deduced from the expression for 
the geometric nonlinearity [7| using the beam thickness 
known from the fabrication and lengths calculated from 
the beam frequencies lo\ , 102 and the material constants 

The response of the system driven near resonance 
and for small dissipation and driving can be calculated 
from ([T]) using the standard methods of secular pertur- 
bation theory [3. This approach has previously used 
for the case of parametrically driven nanomechanical de- 
vices [13, [HI , where hysteretic switches between dif- 
ferent stable states in frequency sweeps were also pre- 
dicted. Briefly, we introduce the slowly varying com- 
plex mode amplitudes Ai^An and forces Fi^Fn using 
X[7] = Re(A[/]e*"[^i*) and .gi5[/] = Re(i^[/] (i)e'"[^i*) (where 
[/] stands for either / or //), substitute into the equa- 
tions of motion, and retain only the near resonant terms. 
This reduces the equations of motion to 

2tu;jAi+iijnjAj+aj\Aj\^Aj+/3i\Aij\^Aj - Fi{t), (2) 

and a corresponding equation for An. The mode nonlin- 



earity parameters a[j] and /3[/] are calculated from ai, a2 
and the eigenvectors e/, e// so that all the parameters in 
(PI are known from linear measurements and the beam 
geometry and material constants. 

We first look at the case where each mode responds at 
the drive frequency ujd[i] which is set near the resonant 
frequency uj[i] so that A[j] oc e*'^'^o[^i~"[^i)*, with 
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[2LUiicjDi - u;i) -ai\Ai\^- l3i\Aii\^f + uj]-fj ' 

(3) 

etc. For a single drive (e.g. An — Fjj — 0), so that 
the cross-mode nonlinear coupling proportional to 
is not involved, this expression reproduces the regular 
Duffing response curve [l6|. Prominent features are the 
shift of the frequency of the maximum response to larger 
values for positive a (nonlinear spring stiffening), and 
bistability and hysteresis that develop above a critical 
drive strength. An experimental example of upward and 
downward frequency sweeps for a drive strength 4.3 times 
critical is shown in Fig. lb. 
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FIG. 2: (color online) Frequency sweeps of the first mode 
response for increasing amplitudes of the second mode: the 
dark (blue) lines are for upward sweeps; the light (yellow) 
line in (a) is a downward sweep, showing the hysteresis and 
amplitude jumps in this case. There are no hysteretic jumps 
for (b) and (c), and so the downward sweeps are not shown. 



We first look at the linear and weakly nonlinear re- 
sponse of the first mode when the second mode is driven 
into its strongly nonlinear regime. Retaining the lead- 
ing order effect of the first mode on the intensity of the 
second mode we can write the numerator in ^ as 



2LoiiujDi - uJi) - I' - ai\Ai\' 



where A^^j is the solution for the second mode in the 
absence of mode I given by the equation for An corre- 
sponding to ([3|) assuming zero Aj. The new parameter 
Q/ is an effective Duffing nonlinearity coefficient and is 
given by a/ = a/ — 2(3il3iiU!ii d\A^j} /duiDii- Equation 
^ predicts two important effects: the frequency tun- 
ing of the first mode proportional to the square of the 
amplitude of the second mode (upwards for positive /?/); 
and the change in the effective nonlinear coefffcient for 
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the motion of the first mode depending on the excitation 
strength of the second mode through the last term in aj. 

To test the frequency tuning we excite the second mode 
at a drive level approximately 4.3 times the critical value 
so that the spectral response is the strongly nonlinear 
Dufhng curve (see Fig. lb). As the actuation frequency 
of the second mode is steadily increased in small steps 
its vibration amplitude rises over a wide frequency range 
until it drops to the lower amplitude state beyond the 
maximum. The evolution of the spectral response of the 
first mode is monitored at a driving level approximately 
four times lower than the critical value for this mode us- 
ing a network analyzer. The dependence of the first mode 
frequency shift on the vibration amplitude of the second 
mode is shown in Fig. Ic. The experimental results for 
frequency tuning closely follow the predicted parabolic 
dependence /?/ I^//]^ for amplitudes up to about 20nm. 
Monitoring the frequency shift of one mode proportional 
to the square of the displacement of a second mode has 
been proposed for quantum nondemolition measurements 
in nanoelectromechanical systems (l7, [l8|. 

If the actuation level of the first mode is increased 
above the onset of nonlinearity, while the second mode 
drive is kept at a much higher level, then the effec- 
tive nonlinear coefhcient a/ is decreased. The decrease 
is largest when the second mode is driven on the por- 
tion of the Duffing response curve where the intensity 
is increasing linearly with the drive frequency |^//|^ — 
2ujii{ujdii — uiii)/aii. This gives the minimum effective 
nonlinear coefficient a/^min = ctj (^1 — ^jfjj^ ■ This re- 
sult indicates that if the coupling is strong enough, so 
that PjPu > ajajj, the minimum value of the effective 
nonlinear coefficient is negative. In this case the reso- 
nance curve tilts to the left, as opposed to the usual case 
for a doubly clamped beam where the peak leans to the 
right. It also means that the nonlinear coefficient van- 
ishes for some drive strength and frequency of the second 
mode. 

Some experimental results for the first mode driven 
at twice the critical strength illustrating this effect are 
shown in Fig. 2. The plots show the shape of the res- 
onance peak for three values of the drive frequency of 
the second mode. In panel (a), the amplitude An of the 
second mode is low and the first mode spectral response 
has the regular nonlinear Duffing shape leaning to the 
right. For larger Ajj as in (b), the first mode resonance 
peak shape assumes a form close to a Lorentzian with 
little nonlinearity apparent. For even larger Ajj as in 
(c) the sign of the effective Duffing coefficient becomes 
slightly negative, causing the spectral response peak to 
lean to the left. The quenching of the nonlinearity in (b) 
could be used to enhance the limited dynamic range ^] 
of nanomechanical devices. 

If the drive level of the first mode is increased further 
a variety of new effects can be observed. Under these 
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FIG. 3: (color online) First mode frequency response as in 
Fig. 2 for increasing values of the second mode drive frequency 
but for stronger driving of the first mode (both modes are 
driven at about four times the critical strength). The top 
plots are experimental measurements while the bottom ones 
are theoretical simulations. The second mode response was 
monitored, but is not shown. 



conditions the dynamics of the system is not fully ex- 
plained by the steady state solutions as in ([3]), and so 
we solve for the expected behavior by numerically inte- 
grating the time-dependent coupled equations Now 
the behavior is more complex, as the response of the first 
mode becomes large enough to cause transitions in the 
second mode response through the nonlinear coupling. 
As a result the spectral response curves acquire peculiar 
nontrivial shapes, as shown in Fig. 3. Numerical sim- 
ulations of the equations give good predictions for the 
complex phenomena. 

An obvious difference between the experimental and 
theoretical plots in Fig. 3 is the noisy regions on the theo- 
retical curves near the up and down transitions. This dif- 
ference actually results from the different ways the plots 
are generated in theory and experiment, and a more care- 
ful investigations shows consistent and interesting dy- 
namics in both experiment and theory: for the drive 
frequencies near the transition points the fast (about 
ITMHz) oscillating response becomes amplitude modu- 
lated with a frequency of about 10 to 20 kHz [i^. The 
theoretical equations ^ show that the frequency for the 
amplitude modulation is determined by the line width, 
which in our experimental setup is about 8 kHz, and is 
not related to sum and difference frequencies of the two 
modes. 

Examples of simulated and measured amplitude mod- 
ulation dynamics for three different sets of drive parame- 
ters are shown in Fig. 4 for a second device with different 
parameters to the one used for Figs. 1-3. The first column 
shows examples of numerically calculated phase portraits 
of Re Aji versus Re Aj obtained by solving the time de- 
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FIG. 4: Complex dynamics of two strongly coupled nanome- 
chanical resonators. The three rows correspond to different 
input parameters (drive frequencies and amplitudes). The 
first column shows the theoretical calculation of the phase 
portrait, the second shows its experimental measurement, and 
the third column shows the corresponding experimental power 
spectrum of the optical measurement near one of the drive 
frequencies. 



pendent equations ([2]). Phase portraits measured experi- 
mentally are shown in the second column |24i] . These are 
obtained using homodyne down-conversion of the trans- 
duced mechanical signals from both modes, which are 
then read by independent oscilloscope channels. We also 
perform wider frequency band spectrum analyzer mea- 
surements shown in the third column. Since the nature 
of the dynamics is sensitive to the precise values of the 
system parameters, the drive strength and frequencies 
in the experiment are slightly adjusted from the values 
corresponding to the theoretical plots to produce compa- 
rable phase portraits. 

The top row of Fig. 4 shows a relatively simple example 
where the motion is periodic but the Aj — An trajectory 
forms a small figure-of-eight loop. The dynamics can be 
roughly understood in terms of transitions between states 
with the first mode at large amplitude and the second 
mode at small amplitude and vice versa. The spectrum 
of the measured mechanical signal shows satellite peaks 
corresponding to the anharmonic amplitude modulation. 

As the parameters of the system are changed, we have 
observed period doubling or quadrupling in the ampli- 
tude modulation, where the phase trajectory takes two 
or four revolutions in order to complete the cycle [l9| . 
An example of period quadrupling is shown in the mid- 
dle row of Fig. 4. The complicated loop structures are 
visible in both theoretical and experimental trajectories, 
and the spectrum reveals amplitude modulation peaks 
at frequencies that correspond to double and quadruple 



periods. Period doubling transitions are often associated 
with chaotic dynamics 19], and indeed for other param- 
eter values we observe chaos in our coupled nanoelec- 
tromechanical system and in the theoretical model, as 
shown in the bottom row of Fig. 4. The evidence for the 
chaotic dynamics in the experiment is the broad band 
component to the spectrum (evident in the shoulders to 
the amplitude modulation peaks), and a phase portrait 
trajectory that does not form a closed loop [25[. The 
theoretical model shows a similar phase portrait. 

Our detailed study of two elastically coupled, indepen- 
dently driven, nanomechanical beam resonators reveals 
complex nonlinear dynamics with a number of potential 
applications. For example, driving one of the modes can 
be used to tune the effective nonlinearity of the other 
mode. This can be used to significantly increase the dy- 
namic range of the resonator by quenching the effective 
nonlinearity. In a first approximation, the motion of one 
mode couples quadratically to the resonance frequency of 
the other mode, a phenomenon that has been proposed 
for quantum nondemolition measurements in nanome- 
chanical systems. For larger vibration amplitudes spon- 
taneous oscillations of amplitude modulation develop, at 
a frequency determined by the resonator ring-down time. 
These oscillations show period doubling and chaos char- 
acteristic of strongly nonlinear systems. The full range 
of complex dynamics investigated is quantitatively re- 
produced by theory. Our success at predicting and sub- 
sequently observing quite delicate features of the non- 
linear dynamics is strong evidence that the nonlinearity 
and coupling in arrays of nanomechanical devices can be 
quantitatively understood and controlled. 

Many of the ideas leading to this work were devel- 
oped in collaboration with Ron Lifshitz supported by the 
U.S. -Israel Binational Science Foundation (BSF) through 
Grant No. 2004339. We thank Matt Mathcny for many 
useful discussions. 
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